Load packages

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.0     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Check working directory

getwd() 
## [1] "/Users/jillashey/Desktop/PutnamLab/Repositories/Astrangia_repo/Astrangia_repo"

Load data

# Load data
daily <- read.csv("data/DailyMeasurements.csv", header = T)
daily <- daily[-13,] # remove tank 14 because it cracked 
daily$Timepoint <- as.Date(daily$Timepoint, "%m/%d/%Y")
daily$Tank <- as.factor(daily$Tank)
daily <- subset(daily, Treatment=="Heat" | Treatment=="Ambient")

Discrete pH calculations from Tris calibrations.

path <-("data/pH_Tris") #set path to calibration files
file.names<-list.files(path = path, pattern = "csv$") #list all the file names in the folder to get only get the csv files
pH.cals <- data.frame(matrix(NA, nrow=length(file.names), ncol=3, dimnames=list(file.names,c("Date", "Intercept", "Slope")))) #generate a 3 column dataframe with specific column names

for(i in 1:length(file.names)) { # for every file in list start at the first and run this following function
  Calib.Data <-read.table(file.path(path,file.names[i]), header=TRUE, sep=",", na.string="NA", as.is=TRUE) #reads in the data files
  file.names[i]
  model <-lm(mVTris ~ TTris, data=Calib.Data) #runs a linear regression of mV as a function of temperature
  coe <- coef(model) #extracts the coeffecients
  summary(model)$r.squared #extracts the r squared
  plot(Calib.Data$mVTris, Calib.Data$TTris) #plots the regression data
  pH.cals[i,2:3] <- coe #inserts coefficients in the dataframe
  pH.cals[i,1] <- substr(file.names[i],1,8) #stores the file name in the Date column
}

colnames(pH.cals) <- c("Calib.Date",  "Intercept",  "Slope") #rename columns
pH.cals #view data
##              Calib.Date  Intercept      Slope
## 20210122.csv   20210122  -86.73848  1.0218902
## 20210126.csv   20210126  -86.23930  0.9322418
## 20210203.csv   20210203  -88.14064  1.0017943
## 20210205.csv   20210205  -88.23311  0.8416322
## 20210217.csv   20210217  -87.53924  1.1237563
## 20210221.csv   20210221  -85.87848  0.8968835
## 20210223.csv   20210223  -86.36221  0.9370372
## 20210225.csv   20210225  -85.86760  0.9234509
## 20210227.csv   20210227  -87.60427  1.0150785
## 20210308.csv   20210308  -89.42691  0.9977500
## 20210317.csv   20210317  -86.21007  0.9799588
## 20210325.csv   20210325  -87.82550  1.0238389
## 20210409.csv   20210409  -88.31599  1.1153461
## 20210411.csv   20210411  -53.62103 -0.1385311
## 20210420.csv   20210420  -68.46091  0.2421390
## 20210524.csv   20210524  -96.95526  1.1877897
## 20210527.csv   20210527  -98.21886  1.2378231
## 20210531.csv   20210531  -88.25076  0.8102227
## 20210603.csv   20210603  -92.83213  1.0612981
## 20210609.csv   20210609  -97.02534  1.2354378
## 20210611.csv   20210611  -97.84557  1.2647915
## 20210615.csv   20210615  -96.14596  1.1781620
## 20210618.csv   20210618  -98.24662  1.2787030
## 20210622.csv   20210622  -96.86842  1.1817890
## 20210625.csv   20210625 -100.50568  1.3608393
## 20210628.csv   20210628  -92.85454  1.0585358
## 20210702.csv   20210702  -95.97090  1.1559734
## 20210706.csv   20210706 -101.52648  1.3059321
## 20210708.csv   20210708  -98.39602  1.1416130
## 20210713.csv   20210713 -102.79400  1.3403794
## 20210720.csv   20210720 -100.01209  1.0809837
#constants for use in pH calculation 
R <- 8.31447215 #gas constant in J mol-1 K-1 
F <- 96485.339924 #Faraday constant in coulombs mol-1

#merge with Seawater chemistry file
SW.chem <- merge(daily, pH.cals, by="Calib.Date")

Calculate total pH.

mvTris <- SW.chem$Temperature*SW.chem$Slope+SW.chem$Intercept #calculate the mV of the tris standard using the temperature mv relationships in the measured standard curves 
STris<-35 #salinity of the Tris
phTris<- (11911.08-18.2499*STris-0.039336*STris^2)*(1/(SW.chem$Temperature+273.15))-366.27059+ 0.53993607*STris+0.00016329*STris^2+(64.52243-0.084041*STris)*log(SW.chem$Temperature+273.15)-0.11149858*(SW.chem$Temperature+273.15) #calculate the pH of the tris (Dickson A. G., Sabine C. L. and Christian J. R., SOP 6a)
SW.chem$pH.Total<-phTris+(mvTris/1000-SW.chem$pH.MV/1000)/(R*(SW.chem$Temperature+273.15)*log(10)/F) #calculate the pH on the total scale (Dickson A. G., Sabine C. L. and Christian J. R., SOP 6a)

Plot daily measurements

# By Treatment
temp.trt = ggplot(SW.chem, aes(x=Treatment, y=Temperature.C)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Temperature°C") +
  theme(axis.text.x = element_text(angle = 90))
temp.trt
## Warning: Removed 11 rows containing non-finite values (stat_boxplot).

sal.trt = ggplot(SW.chem, aes(x=Treatment, y=Salinity.psu)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Salinity.psu") +
    theme(axis.text.x = element_text(angle = 90))
sal.trt
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).

pH.trt = ggplot(SW.chem, aes(x=Treatment, y=pH.Total)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("pH Total Scale") +
  theme(axis.text.x = element_text(angle = 90))
pH.trt
## Warning: Removed 24 rows containing non-finite values (stat_boxplot).

light.trt = ggplot(SW.chem, aes(x=Treatment, y=Light)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Light µmol photons m^2 s^-1") +
    theme(axis.text.x = element_text(angle = 90))
light.trt
## Warning: Removed 1059 rows containing non-finite values (stat_boxplot).

flow.trt = ggplot(SW.chem, aes(x=Treatment, y=Flow.mL.sec)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Flow mL/sec") +
    theme(axis.text.x = element_text(angle = 90))
flow.trt
## Warning: Removed 1760 rows containing non-finite values (stat_boxplot).

plot.trt <- grid.arrange(temp.trt, sal.trt, pH.trt, light.trt, flow.trt, ncol=3, nrow = 2, clip="off")
## Warning: Removed 11 rows containing non-finite values (stat_boxplot).
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).
## Warning: Removed 24 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1059 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1760 rows containing non-finite values (stat_boxplot).

ggsave("output/Daily_Measurements.By.Treatment.pdf", plot.trt, width = 21, height = 21, units = c("in"))

# By Tank
temp.tank = ggplot(SW.chem, aes(x=Tank, y=Temperature.C)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Temperature°C") +
    theme(axis.text.x = element_text(angle = 90))
temp.tank
## Warning: Removed 11 rows containing non-finite values (stat_boxplot).

sal.tank = ggplot(SW.chem, aes(x=Tank, y=Salinity.psu)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Salinity.psu") +
    theme(axis.text.x = element_text(angle = 90))
sal.tank
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).

pH.tank = ggplot(SW.chem, aes(x=Tank, y=pH.Total)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("pH Total Scale") +
    theme(axis.text.x = element_text(angle = 90))
pH.tank
## Warning: Removed 24 rows containing non-finite values (stat_boxplot).

light.tank = ggplot(SW.chem, aes(x=Tank, y=Light)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Light µmol photons m^2 s^-1") +
    theme(axis.text.x = element_text(angle = 90))
light.tank
## Warning: Removed 1059 rows containing non-finite values (stat_boxplot).

flow.tank = ggplot(SW.chem, aes(x=Tank, y=Flow.mL.sec)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Flow mL/sec") +
    theme(axis.text.x = element_text(angle = 90))
flow.tank
## Warning: Removed 1760 rows containing non-finite values (stat_boxplot).

plot.tank <- grid.arrange(temp.tank, sal.tank, pH.tank, light.tank, flow.tank, ncol=3, nrow = 2, clip="off")
## Warning: Removed 11 rows containing non-finite values (stat_boxplot).
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).
## Warning: Removed 24 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1059 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1760 rows containing non-finite values (stat_boxplot).

ggsave("output/Daily_Measurements.By.Tank.pdf", plot.tank, width = 21, height = 21, units = c("in"))

test for differences between tanks and treatments

# By Treatment
mod.temp <- aov(Temperature.C ~ Treatment, data = SW.chem)
mod.temp
## Call:
##    aov(formula = Temperature.C ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares    6140.63  53132.89
## Deg. of Freedom         1      1827
## 
## Residual standard error: 5.392777
## Estimated effects may be unbalanced
## 11 observations deleted due to missingness
summary(mod.temp)
##               Df Sum Sq Mean Sq F value Pr(>F)    
## Treatment      1   6141    6141   211.1 <2e-16 ***
## Residuals   1827  53133      29                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 11 observations deleted due to missingness
hist(mod.temp$residuals)

mod.sal <- aov(Salinity.psu ~ Treatment, data = SW.chem)
mod.sal
## Call:
##    aov(formula = Salinity.psu ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares    17.6973 2115.0465
## Deg. of Freedom         1      1815
## 
## Residual standard error: 1.079498
## Estimated effects may be unbalanced
## 23 observations deleted due to missingness
summary(mod.sal)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Treatment      1   17.7  17.697   15.19 0.000101 ***
## Residuals   1815 2115.0   1.165                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 23 observations deleted due to missingness
hist(mod.sal$residuals)

mod.pH <- aov(pH.Total ~ Treatment, data = SW.chem)
mod.pH
## Call:
##    aov(formula = pH.Total ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares    4.32348  91.46740
## Deg. of Freedom         1      1814
## 
## Residual standard error: 0.2245508
## Estimated effects may be unbalanced
## 24 observations deleted due to missingness
summary(mod.pH)
##               Df Sum Sq Mean Sq F value Pr(>F)    
## Treatment      1   4.32   4.323   85.74 <2e-16 ***
## Residuals   1814  91.47   0.050                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 24 observations deleted due to missingness
hist(mod.pH$residuals)

mod.light <- aov(Light ~ Treatment, data = SW.chem)
mod.light
## Call:
##    aov(formula = Light ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares        0.6  380831.6
## Deg. of Freedom         1       779
## 
## Residual standard error: 22.11046
## Estimated effects may be unbalanced
## 1059 observations deleted due to missingness
summary(mod.light)
##              Df Sum Sq Mean Sq F value Pr(>F)
## Treatment     1      1     0.6   0.001  0.972
## Residuals   779 380832   488.9               
## 1059 observations deleted due to missingness
hist(mod.light$residuals)

mod.flow <- aov(Flow.mL.sec ~ Treatment, data = SW.chem)
mod.flow
## Call:
##    aov(formula = Flow.mL.sec ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares     3167.9  679148.5
## Deg. of Freedom         1        78
## 
## Residual standard error: 93.31148
## Estimated effects may be unbalanced
## 1760 observations deleted due to missingness
summary(mod.flow)
##             Df Sum Sq Mean Sq F value Pr(>F)
## Treatment    1   3168    3168   0.364  0.548
## Residuals   78 679148    8707               
## 1760 observations deleted due to missingness
hist(mod.flow$residuals)

# By Tank
mod.temp <- aov(Temperature.C ~ Tank, data = SW.chem)
mod.temp
## Call:
##    aov(formula = Temperature.C ~ Tank, data = SW.chem)
## 
## Terms:
##                     Tank Residuals
## Sum of Squares   6154.91  53118.61
## Deg. of Freedom       15      1813
## 
## Residual standard error: 5.412831
## Estimated effects may be unbalanced
## 11 observations deleted due to missingness
summary(mod.temp)
##               Df Sum Sq Mean Sq F value Pr(>F)    
## Tank          15   6155   410.3   14.01 <2e-16 ***
## Residuals   1813  53119    29.3                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 11 observations deleted due to missingness
hist(mod.temp$residuals)

mod.sal <- aov(Salinity.psu ~ Tank, data = SW.chem)
mod.sal
## Call:
##    aov(formula = Salinity.psu ~ Tank, data = SW.chem)
## 
## Terms:
##                      Tank Residuals
## Sum of Squares    18.8038 2113.9400
## Deg. of Freedom        15      1801
## 
## Residual standard error: 1.083402
## Estimated effects may be unbalanced
## 23 observations deleted due to missingness
summary(mod.sal)
##               Df Sum Sq Mean Sq F value Pr(>F)
## Tank          15   18.8   1.254   1.068  0.382
## Residuals   1801 2113.9   1.174               
## 23 observations deleted due to missingness
hist(mod.sal$residuals)

mod.pH <- aov(pH.Total ~ Tank, data = SW.chem)
mod.pH
## Call:
##    aov(formula = pH.Total ~ Tank, data = SW.chem)
## 
## Terms:
##                     Tank Residuals
## Sum of Squares   4.57224  91.21864
## Deg. of Freedom       15      1800
## 
## Residual standard error: 0.2251156
## Estimated effects may be unbalanced
## 24 observations deleted due to missingness
summary(mod.pH)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Tank          15   4.57 0.30482   6.015 2.01e-12 ***
## Residuals   1800  91.22 0.05068                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 24 observations deleted due to missingness
hist(mod.pH$residuals)

mod.light <- aov(Light ~ Tank, data = SW.chem)
mod.light
## Call:
##    aov(formula = Light ~ Tank, data = SW.chem)
## 
## Terms:
##                     Tank Residuals
## Sum of Squares   17193.7  363638.5
## Deg. of Freedom       15       765
## 
## Residual standard error: 21.80239
## Estimated effects may be unbalanced
## 1059 observations deleted due to missingness
summary(mod.light)
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## Tank         15  17194  1146.2   2.411 0.00199 **
## Residuals   765 363638   475.3                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1059 observations deleted due to missingness
hist(mod.light$residuals)

mod.flow <- aov(Flow.mL.sec ~ Tank, data = SW.chem)
mod.flow
## Call:
##    aov(formula = Flow.mL.sec ~ Tank, data = SW.chem)
## 
## Terms:
##                     Tank Residuals
## Sum of Squares   51655.4  630661.0
## Deg. of Freedom       15        64
## 
## Residual standard error: 99.26771
## Estimated effects may be unbalanced
## 1760 observations deleted due to missingness
summary(mod.flow)
##             Df Sum Sq Mean Sq F value Pr(>F)
## Tank        15  51655    3444   0.349  0.987
## Residuals   64 630661    9854               
## 1760 observations deleted due to missingness
hist(mod.flow$residuals)

Plot by date

# Plot and save
# By Treatment
temp.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=Temperature.C, group=Treatment)) +
  geom_line(aes(color = Treatment), size = 1) 
sal.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=Salinity.psu, group=Treatment)) +
  geom_line(aes(color = Treatment), size = 1) 
ph.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=pH.Total, group=Treatment)) +
  geom_line(aes(color = Treatment), size = 1) 
light.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=Light, group=Treatment)) +
  geom_line(aes(color = Treatment), size = 1) 
flow.trt.date <- ggplot(SW.chem, aes(x=Timepoint, y=Flow.mL.sec, group=Treatment)) +
  geom_line(aes(color = Treatment), size = 1) 

plot.trt.timepoint <- grid.arrange(temp.trt.date, sal.trt.date, ph.trt.date, light.trt.date, flow.trt.date, ncol=1, nrow = 5, clip="off")
## Warning: Removed 16 row(s) containing missing values (geom_path).
## Warning: Removed 1168 row(s) containing missing values (geom_path).

ggsave("output/Daily_Measurements.By.Date.Treatment.pdf", plot.trt.timepoint, width = 30, height = 21, units = c("in"))

# By Tank
temp.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=Temperature.C, group=Tank)) +
  geom_line(aes(color = Tank), size = 1) 
sal.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=Salinity.psu, group=Tank)) +
  geom_line(aes(color = Tank), size = 1) 
ph.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=pH.Total, group=Tank)) +
  geom_line(aes(color = Tank), size = 1) 
light.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=Light, group=Tank)) +
  geom_line(aes(color = Tank), size = 1) 
flow.tank.date <- ggplot(SW.chem, aes(x=Timepoint, y=Flow.mL.sec, group=Tank)) +
  geom_line(aes(color = Tank), size = 1) 

plot.tank.timepoint <- grid.arrange(temp.tank.date, sal.tank.date, ph.tank.date, light.tank.date, flow.tank.date, ncol=1, nrow = 5, clip="off")
## Warning: Removed 16 row(s) containing missing values (geom_path).

## Warning: Removed 1168 row(s) containing missing values (geom_path).

ggsave("output/Daily_Measurements.By.Date.Tank.pdf", plot.tank.timepoint, width = 21, height = 25, units = c("in"))

plot.all <- grid.arrange(temp.trt.date, temp.tank.date, sal.trt.date, sal.tank.date, ph.trt.date,ph.tank.date,  light.trt.date, light.tank.date, flow.trt.date, flow.tank.date, ncol=2, nrow=5)
## Warning: Removed 16 row(s) containing missing values (geom_path).
## Warning: Removed 16 row(s) containing missing values (geom_path).
## Warning: Removed 1168 row(s) containing missing values (geom_path).

## Warning: Removed 1168 row(s) containing missing values (geom_path).

ggsave("output/Daily_Measurements.By.Date.pdf", plot.all, width = 30, height = 25, units = c("in"))